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ABSTRACT 

We present a new hybrid Smoothed Particle Hydrodynamics (SPH)/7V-body method 
for modelling the collisional stellar dynamics of young clusters in a live gas background. 
By deriving the equations of motion from Lagrangian mechanics we obtain a formally 
conservative combined SPH- TV-body scheme. The SPH gas particles are integrated 
with a 2nd order Leapfrog, and the stars with a 4th order Hermite scheme. Our new 
approach is intended to bridge the divide between the detailed, but expensive, full hy- 
drodynamical simulations of star formation, and pure iV-body simulations of gas-frcc 
st ar clusters. We ha ve implemented this hybrid approach in the SPH code SEREN 



(jHubber et al.ll2011l ) and perform a series of simple tests to demonstrate the fidelity 
of the algorithm and its conservation properties. We investigate and present resolu- 
tion criteria to adequately resolve the density field and to prevent strong numerical 
scattering effects. Future developments will include a more sophisticated treatment of 
binaries. 
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1 INTRODUCTION 

The formation and dissolution of young stellar clusters is an important, but complex problem that requires computer simula- 
tions to_exrjlc£e ind^aiL_Stai^ of tens to thousands of stars 



(e.g. lElmegreenl I2OO0I : lLada fc Ladall2003l : iMcKee fc Ostrikerl 120071 : iKlessen et aPboogl ). Therefore the early phases of the 



dynamical evolution of star clusters involves young stars moving through a significant and dynamic gaseous background. To 
understand the complete process of cluster formation and evolution requires a combination of self-gravitating hydrodynamics 
for the gas from which stars and planets form, and the gravitational iV-body dynamics of (multiple) stars and planets once 
they have formed. 

Previously, detailed hydrodynamics and accurate iV-body dynamics have been separated. Hydrodynamical simulations 
have been used to simulate the turbulent gas dynamics leading to fragmentation and star formation, while iV-body simulations 
tend to follow the late gas-free stages of star cluster evolution. However, there is a very signifcant and important phase in the 
life of a star cluster in which stellar dynamics within a gas background is vitally important. In the 'gas-rich' phase, which 
occurs around 1-5 MyiQ the stars are interacting dynamically in a live gas background. The star formation process tends to 
produce binary and multiple systems in complex hiearchical structures. Dynamical interactions between single and multiple 
systems during the su bsequent few Myr changes the bina ry properties of the stars as well as the structure and dynamics of 
the whole cluster (see Allison et al. 20091 ; Goodwin 2010l . and references therein). Therefore, the binary pro perties of star s 
released into the field after gas expulsion will depend on stellar dynamics during the gas-rich phase (see also lKroupalll995l ). 
In addition, the early stages of planet formation will occur during this g as-rich phase and interactions may seriously alter the 
architecture and properties of planetary systems I Parker fc Quan3201ll ). Accurate observations of the binary and dynamical 



1 We note that star formation is not an instantaneous process and that stars are still forming whilst others are interacting and dynamically 
evolving. However, we make a rough first approximation that most stars form in the first Myr, they then evolve in a dynamical gas 
background which is expelled at 5 Myr. 
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properties of clusters are usually only available once gas is expelled (especially those to be provided by Gaia), which means 
they will have been altered by dynamical evolution in the gas-rich phase. 

In hydrody namical simulations, we genera lly replace the dense collapse pha se of gas into stars wi th sink particles (see 
Bate et alii 19951 for SPH; Krumholz et al. 2004 for AMR impl ementations, see also Federrath et al.ll201C ). These sink particles 



can r epresent individual stars if their sizes are < 1 AU (e.g. bate et al.ll2003l : TGoodwin et al.ll2004l : lBatell2009l : Offner et al 



20091') or larger regions perhaps containing primordial mult iple systems which cannot be resolved if their sizes are > 10 AU (e.g 



Bonnell et akll2004l : ISmith et al1l2009l : Ijappsen et al]|2005l ). Sinks have the huge advantage of allowing dense, computationally 



expensive, and unresolvable regions to be 'compressed' into a particle which can interact with the surrounding gas and accrete 
from it. However, sink particles are not point-like A r -body particles as, even if each sink represents a single star, (a) their 

gravity is softened, and (b) sinks accrete from the surrounding gas. 

Pure ./V-body simulations of stellar systems have a long history fsee lAarsethlbooj iHeggie fe Hutll2003h . but most ignore 
the early gas-rich phases of a star cluster's life. The usual way to incl ude gas and model the gas-rich phase is to introduce an 
external potential, which is often a simple Plumm er or King mo del (e.g. Lada et al.ll984 ; Goodwin1ll997 ; Baumgardt fe Kroupa 
2007 ; Moeckel fe Clarke] 2011 ; Smith et al. 2011 ). although see Gever fc Burkert (|2001 ) used softened 'star particles' in a live 



gas background. Generally, the external potential is allowed to vary with time, e.g. to model the expulsion of gas from a 
cluster. However, the use of a simple analytic external potential to model the gas (which is often the majority of the mass in 
an embedded cluster) is clearly a vast over-simplification. 

In this paper we introduce a new hybrid iV-body/Sm oothed Particle Hydrodynamics (SPH) algorithm that has been 
implemented in the SPH code SEREN (jHubber et al.ll201ll ). The stellar dynamics are computed with a 4th-order integrator 
allowing the details of iV-body interactions between stars to be followed. SPH gas particles are used to represent a live 
background gas potential in which the stellar dynamics is modelled. We emphasise that this hybrid method is not a replacement 
for fully self-consistent, high-resolution star formation simulations or detailed iV-body simulations. Rather, it represents a fast 
way of exploring stellar dynamics in a live background potential which can be used to perform large suites of simulations to 
explore large parameter space, or to inform the initial conditions of pure iV-body simulations of the post-gas phase. 

In Section O we introduce the hydrodynamical and N-body methods used and how they are combined algorithmically. 
In Section [3] we present a number of simple tests to demonstrate the accuracy and robustness of our method. In Section [4] 
we discuss various important caveats of our method, in particular understanding resolution effects, and also discuss possible 
astrophysical problems that can be explored with this code. 



2 NUMERICAL METHOD 



Self-gravitating hydrodynamical simulations in astrop hysics are u s ually modelled using either a Lagrangian, particle-based 
approach such as Smoothed Particle Hydrodynamics l|Lucvlll977 ; Gingold fc Monaghanl Il977l ). or an Eulerian, grid-based 
approach such as Adaptive Mesh Refinement Hydrodynamics ( Berger fc ColellJ 1989h ■ Whereas SPH derives interaction 
terms by computing particle-particle force terms and integrating the motion of each particle individually, grid codes operate 
by computing fluxes across neighbouring grid cells. Since N-body codes also work by computing forces and integrating positions 
and velocities, SPH is the most natural hydrodynamical method to merge directly with TV-body dynamics as the particle- 
nature of the gas and stars are easily compatible making it straight-forward to derive the coupling force terms and to merge 
their individual integration schemes. 



We use a conservative self-gravitating SPH formulation (|Price fc Monaghanll2007l ) to model the gas dynamics and include 
the star particles within the SPH formulation as a s pecial type of SPH particle , rather than external iV -body particles. 
Following most modern conservative SPH schemes (e.g. Springel fc Hernquist|[20o3 ; Price fc Monaghanll2007l ). the smoothing 
length of a gas particle i is set by the relation, 

* = ■ <» 

and the SPH gas density is given by 

N 



(2) 



where r^, hi, rrii, pi are the position, smoothing length, mass and density of particle i respectively, 



r,, W is the 



SPH smoothing kernel and r\ is a dimensionless number that controls the mean number of neighbours (usually set to 1.2 
to have ~ 60 neighbours). Since hi and pi depend on each other, we must iterate between Equations Q] and [2] in order to 
reach a consistent solution. In contrast, the star particles have a constant smoothing length which represents the gravitational 
softening length to p revent violent 2-body collisions with other stars, in place of using more complicated algorithms such 
as regularisation (See Aarseth 20031 . for a description of common N-body techniques). In order to reduce 'scattering' during 
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star-gas interactions, we use the mean-smoothing length approach (See Appendix A of iPrice fc Monaghanl 120071 ). to keep 
star-gas interactions as smooth as possible. We can now formulate the Lagrangian of the system containing all interaction 
terms and then derive the equations of motion via the Euler-Lagrange Equations. This simple approach allows us to develop 
a conservative scheme which in principle can be integrated to arbitrary accuracy (i.e. if direct summation of gravitational 
forces and a constant, global timestep is used) . Due to the larger energ y errors often produced by N-body encounters, we 
use a higher-order Hermite integration scheme ( Makino fc Aarsethl 19921 ) to integrate star particles, and a simpler 2nd-order 
Leapfrog kick-drift-kick scheme to integrate the gas particles motion. 



2.1 Gravitational force softening in SPH 



The gravitati onal force softening between SPH particles can be derive d in a number of ways (e.g. Plummer softening, See 
DehnenllioOll ). However, it has been suggested bv lBate fc Burkertl (|1997l ) that it is safest to use the SPH kernel itself to derive 



the softening terms to prevent artificial gravitational fragmentation. They showed that for gas condensations where the Jeans 
length was of order the smoothing length or greater than, the net hydrodynamical force is stronger than the net gravitational 
force from all neighbouring particles, thereby suppressing or even reversing the collapse of the condensation and preventing 
fragmentation. We theref ore derive the softening terms from the SPH kernel following the method and nomenclature of 
Price fc Monaghanl (|2007h . 

First, we consider the case of uniform smoothing length. The gravitational potential at the position of particle i due to 
a distribution of SPH particles is given by 



m b i 



(r ab , h) , 



(3) 



where <f> is the gravitational softening kernel (jPrice fc Monaghanl 120071 ) and h is the smoothing length of all SPH particles. 
We note that Eqn. [3] requires the softening kernel to be a negative quantity. The potential is related to the density field by 
Poisson's Equation, 



V 2 $(r) 



4irGp(r) 



(4) 



where the SPH density defined at r a is given by Eqn. [2] This allows us to directly relate the softening kernel to the SPH 
smoothing kernel for a consistent formulation of self-gravity in SPH. Substituting Equations [2] & [3] into Poisson's Equation, 
we obtain 



W{r,h) = 



47rr 2 dr 



(5) 



By direct integration of Eqn.[5]with the appropriate limits, we can obtain the gravitational softening kernel via the gravitational 
force kernel, cj>' where 



<t>'(r, h) = ^(r, h) = % I W(r', h) r' 2 dr' 
or r z 



(6) 



and 



/ r r /-en > 

<f>{r,h) = 47r I ^ W(r',h)r' 2 dr' + J W (r , h) r dr' - J W(r',h)r'dr' 



(7) 



where 1Z is the compact support of the kernel (e.g. 1Z — 2 for the M4-kernel). 

When using variable smoothing lengths, we have two choices for symmetrising the gravitational interaction; (a) use the 
average of the two softening kernels, or (b) use the mean smoothing length in the softening kernel, i.e. 



* a = G E mb ' 

f>_i 



{r ab , h a ) + (p(r ab , h b ) 



or <jy = G ^ m b <j>(r ab , h ab 



(8) 



where h ab = | (h a + h b ). IPrice fc Monaghanl (|2007l ) advocate using the average softening kernel approach since it requires 



less loops over all particles than the mean smoothing length approach. When using only SPH particles where the smoothing 
length is determined by Eqn. [T] there is little difference in the results since the two methods give similar potentials and 
forces. However, when we include star particles which can have an arbitrary small smoothing length, there can be very large 
discrepancies between the two methods. For example, consider a 'collision' between a gas particle and an SPH particle, i.e. 
where the two particle lie at almost the same position in space. The mean softening kernel approach has two terms, one 
which will be quite small due to the smoothing length of the gas particle, and the second term using the star smoothing 
length which can become very large producing a corresponding large scattering force. The mean smoothing length approach 
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however, can never produce a large scattering force, even if the smoothing length of the star becomes zero since the average 
of the two smoothing lengths can never be less than | hi or \ hj . Furthermore, the mean smoothing length method allows us 
to use star particles with zero smoothing length, permitting the study of truly collisional stellar dynamics with unsoftened 
star-star forces, but softened star-gas interactions. Therefore, for the case of our hybrid formulation, we advocate using the 
mean smoothing length approach and we derive all subsequent equations of motion using this method. 



2.2 Conservative SPH 



The SPH equations of motion can be derived from Lagrangian mechanics, resulting in a set of equations that automati- 
cally conserve linear momentum, to machine precisiom and angular mom e ntum and energy, both to integration error (See 
Springel fc Hernquistl [20021 ; iPrice fc Monaghanl 120071 ). iPrice fe Monaghanl (|2007h derived the equations of motion for self- 
gravitating systems with variable smoothing lengths using Lagrangian mechanics. Following their method, we derive the 
equations of motion for a set of SPH particles (with variable smoothing length given by Eqns. [1] & [2j and stars (with a fixed 
smoothing length). If we have N g gas particles with labels b — 1,2, N g and N 3 star particles with labels i = 1,2, N s , 
then by inserting all terms into the Lagrangian, we obtain 



Ni 



C — - ^2 m b vl + - m i v 'i — mb Ub + ^ 



(9) 



where C, 



is the gravitational contribution to the Lagrangian given by 



G 



m b m c 4> bc (h bc ) - ^-^2 m i m j<t>ij(h 



Ng N 3 

G^2 m b m i ( Pbi(h b i) 



(10) 



6=1 c=l i — 1 j — 1 6=1 i = l 

We note that throughout this paper, summations over all SPH gas particles are given by the indices b, c and d and summations 
over all star particles by i, j and k. The equations of motion can be obtained by inserting this Lagrangian into the Euler- 
Lagrange equations, 

d_ fdC \ _ dC_ 
dt \ dv a J dr a 

After inserting the correct terms and evaluating the algebra (See Appendix [A] for a full derivation), we obtain the following 
equation of motion for SPH gas particles 

a fJ = - ^2 m b 



6=1 



0. 



Pa dW ab ,, , , P b dW ab ., 

op a„ (V 



(11) 



p 2 a Q, a dr a 



pffib dr a 



N, 



m b 



a, 



dr a 



(ha) 



-G"^2 m b 4>'ab{hab) r ab -G^rrii (t>'ai{h a i)v ai - — ^ 
6=1 i=l 6=1 

where P a = (7 — 1) p a u a is the thermal pressure of particle a, and Q a , (a, and \ a are defined by 



(Co + Xa) dW ab 



(Cb + Xb) OWal 



dv a 



(In, 



„ , dh a \ -» dW ab 
U a = 1 - - — y m b (ha) . 

opa f-[ ah 



(12) 



(13) 



J dh a s-^ 84>ab ,-r s 

dp a f-1 dh ab 



(14) 



dha °4>ai f-r 

Xa = -t— ym iT =^(hai 

dpa frt Ohai 



(15) 



The Q term is the familiar 'grad-h' correction term that appears in conservative SPH with vary ing smoothing length (e.g . 
Springel fc Hernquist|[2002 ; Price fc Monaghan 2007 ). The £ term is the correction term derived by Price fc Monaghan (2007) 
for gravitational interactions between gas particles in conservative SPH. We obtain an analogous correction term for the 
star-gas interaction, x, which is a summation over all neighbouring star particles. However, \ i s s l m included in a summation 
over all neighbouring gas particles since it is the variation in the smoothing length (which is determined by neighbouring 
particle positions) that gives rise to the correction terms. We have some choice in how to evolve the thermal properties of the 
gas particles. We chose to evolve the specific internal energy equation, i.e. 



du a 

~dt~~pjn 



Pa ST- dWal 



(16) 
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For the star particles, we obtain the following equations of motion for star s, 



N„ 



&s = -G^mi (j>' st {hsi) r si - G^ m b (j>' sb (h sb ) r sb . (17) 



6=1 



We note that since only the smoothing lengths of the gas particles are allowed to vary, all correction terms derived via the 
Lagrangian appear in the equations of motion for the gas particles. 



2.3 Coupled-integration scheme 



The equations of motion for both gas and star particles can be integrated either with a single integration scheme, or with 
two independent integration schemes in parallel. Current SPH codes typically use 2nd-order schemes, such as the Leapfrog 
or the Runge-Kutta-Fehlberg, whereas N-body codes use at least 4th-ord er schemes suc h as the Hermite scheme. In our 
implementation, we chose to use a 2nd-order Leapfrog kick-drift-kick schem e ]Springelll2005l ) to integrate the SPH gas particles 
coupled with a 4th-order Hermite integration scheme (jMakino fc Aarsetblll992l ) to integrate the star particles. One important 
reason for t his choice is that a 4th-order Hermite scheme can be considered as the higher-order equivalent of the leapfrog 
scheme (See lHut et al.l ll995l). where the force, prediction and correction steps are all computed at the same points in the 
timestep for both schemes. We discuss the details of our implementation of both integration schemes, and in particular we 
discuss the modifications to the 4th-order Hermite scheme to include SPH smoothing. 



2.3.1 2nd-order SPH leapfrog integration scheme 

We integrate the SPH particles using a 2nd-order Leapfrog kick-drift-kick integration scheme. A traditional leapfrog works 
by advancing the positions and velocities half-a-step apart, i.e. 

< +1/2 = vr 1/2 +a^At, (18) 
r n a +1 = r n a +v: +1/2 At. (19) 

It is possible to transform the traditional leapfrog equations into a form where the positions and velocities are both updated 
at the end of the step, 

v n a +1 = r^+v^At + ia^At 2 , (20) 

n+1 n . 1 / n , n-\~l\ \ , /r»i\ 

v a = v a + - (a a + a a ^ ) At . (21) 

This form of the leapfrog (also known as the Velocity- Verlet integration scheme) has 3rd-order accuracy in integrating the 
positions and 2nd-order accuracy in integrating the velocities, and therefore 2nd-order overall. As we will see in Section [2.3.21 
it also has the useful property of having its acceleration and 'correction' terms calculated in sync with the corresponding 
terms for the 4th order Hermite scheme used for integrating the star particles. Other integration schemes (e.g. the 2nd-order 
Leapfrog drift-kick-drift) do not necessarily share this property. 

T he timesteps are dete rmined by taking the minimum of three separate conditions, the SPH Courant-Friedri chs-Lewy con 



dition ( Courant et al.l[l928h . the acceleration condition, and the energy condition when cooling is employed (See Hubber et al 



2011). 



7cOUR feg ( 7ACCEL hq \ TeWERGY M ° 1 (r,r)\ 

Ml + 1.2a) c a + (1 + 1.2/3) fc |V-v a | ' ^ |a a |+ e J ' \u - a | + e J ^ 

where the 7 terms are dimensionless timestep multipliers and e ^ 1 is a small number used to prevent a divide-by-zero in the 
case of laJ =0 or \u a \ = 0. 



2.3.2 ^th-order N-body Hermite integration scheme 



We use a standard 4th-ord er Hermite integrat or (|Makino fc Aarsethlll992f ) modified by including the same SPH softening 
scheme as the gas particles (jHubber et al.ll201lT ) to integrate the motion of the star particles. In the Hermite scheme, the first 
time derivative of the acceleration (often referred to as the 'jerk') must be calculated explicitly from Equation [TT] By taking 
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Force, Monopole --0- - 
Force, Quadrupole - -•- - 
Jerk, monopole — v — 




0.01 0.1 



max 

Figure 1. The RMS fractional error using tree gravity of (a) gravitional force with monopole terms only (blue diamonds), (b) gravitational 
force including the quadrupole terms (solid black circles), and (c) gravitational 'jerk' with monopole terms only (red triangles) using the 
geometric MAC as a function of opening angle, $ MAX . 



the time derivative of Equation [17] and using Equation [5l we obtain for the jerk term, 



en Vs! + 3G /L rsi ~ 4nG l^ 



mi (r si ■ v ai ) W{v 3 i, h ai ) 



i=l 
N 



i = l 
N 



•s-^ m b (j)'(T sb ,h s b) , „ „ m b ( r sb ■ v sb ) 0'(r 36 , h ab ) 

■ Or i i V sb + O Or 1 T 3b 



i=l 
N 



r ab 



m b (r ab ■ v s b) W{r 3i , h ab ) 



r ab . (23) 

6=1 '"-I b=1 i— i b=1 l^fl 

Since the jerk can be computed from a single sum over all particles, we can compute it explicitly at the same time as computing 
the accelerations. Once both terms are computed, we can calculate the predicted positions and velocities of the stars at the 
end of the steps, i.e. 



Tl+l n . n A _r i TO a .2 . 1 - TO . ,3 

r, = r a + v a At + -a a At + -a a At . 



6 



n + l 



71 71 a i 1.T7A.2 

r, + a e At H — a. At . 



(24) 
(25) 



The acceleration and jerk are then recomputed at the end of the step, a" +1 and a n+1 . This then allows us to compute the 
higher order time derivatives, 



-3(a£ 



(2a? + a? +1 )At) 



At 2 



6(2(a?-a?+ 1 ) + (a?+a? +1 )At) 
Finally, we apply the correction step where the higher-order terms are added to the position and velocity vectors, i.e. 



71 + 1 



, -U? At 4 + — a™ At 5 
24 3 120 

+ ia n At 3 + — a? At 4 . 
6 24 



We compute the N-body timesteps using the Aarseth timestep criterion 1 Aarseth 20031 ) . 
At a = 7 S 



a* a J + a s 



(26) 
(27) 

(28) 
(29) 

(30) 



where 7 S i s the timestep multiplier for stars. For the very first timestep, we must compute the 2nd and 3rd time derivatives 
explicitly 1 Aarseth 20031 ) since we do not yet have information on the 2nd and 3rd derivatives (since they are only first 
computed at the end of the first timestep). Hereafter, we use Eons. [26l fc 1271 for computing these derivatives. 



2.4 Calculating gravitational terms 

SEREN ( Hubber et al. 2011 ) uses a Barnes-Hut gravity tree ( Barnes fc Hut 19861 ) for computing gravitational forces for all 
self-gravitating gas particles. We use the same tree for computing the gravitational forces due to all gas particles for both 
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SPH gas particles, and star particles. SEREN has a variety of tree-opening criteria that c an be selected at co mpilation-time. 
For most simulations in this paper, we use the Eigenvalue multipole-acceptance criterion (jHubber et al.ll201lT ) because it has 
better force error control, and therefore ultimately better energy error control, than the standard geometric opening-angle 
criterion often used. For nearby SPH particles, which require direct computation of the gravitational acceleration, we also 
compute the jerk term when computing the acceleration of star particles. For tree cells, we compute the jerk contribution due 
to the centre of mass of the cell; however, we ignore the contribution due to the quadrupole moment terms for simplicity. For 
forces due to star particles, we sum all the contributions for the gravitational acceleration (and jerk) directly without using 
a gravity tree. 

In Figure [1] we plot the RMS force (monopole and quadrupole) and jerk (monopole only) errors using the geometric 
opening-angle criterion (in order to plot both monopole and quadrupole errors) for stars in a star-gas Plummer sphere (as 
discussed in Section f3.2[) . We can see that the jerk error scales in the same way as the monopole-only force error. The force 
error using quadrupole moments scales much better than both monopole errors, as expected. In order to allow high-accuracy in 
the calculation of the jerk while still using the tree, we use two different tree opening criteria; one for SPH gas particles (which 
do not need to calculate the jerk), and a stricter one for N-body particles that use the tree. It should be noted that Figure 
[1] represents an upper limit to the expected jerk error. The dominant contribution to the jerk will be from close encounters 
with other stars, which is computed exactly. 



2.5 Block timestepping 

SEREN l|Hubber et alJboilh uses a standard block-timestepping scheme used in many N-body and SPH codes. The timesteps 
of all gas and star particles are restricted to being At = At MAX / 2 n where n — 0, 1, 2, 3, 4, .. is a positive integer. All particles 
and stars therefore are all exactly synchronised on the longest timestep, when the timestep level structure is recomputed. In the 
standard block-timestepping scheme, particle timesteps are only computed once their current timestep has been completed. 
At the end of the step, particles are allowed to move to any lower timestep (higher n), or are allowed to move up one level 
(lower n) provided the new higher level is synchronised with the old level. This approach means particles can rapidly and 
immediately drop to short timesteps when required, but are only allowed to rise back up slowly to prevent timesteps oscillating 

up and down too frequently. 

SEREN also contains the neighbour-timestep monitoring procedure of Saitoh fc Makino ( 2009^ to prevent large timestep 
differences between SPH neighbours generating large energy, momentum and angular momentum errors. Although this algo- 
rithm is not necessarily needed for the tests presented in this paper, it will almost certainly be required for future applications 
where feedback processes can suddenly generate large discontinuities in density and temperature, resulting in large timestep 
disparities. 



3 TESTS 

We present a small suite of numerical tests demonstrating the ac curacy of our hybrid SPH/iV-body approach. Tests using the 



SPH and JV-body components independently were presented bv lHubber et al.l (|2011h . The SPH component was tested using 
a variety of shock-tube, Kelvin-Helmholtz instability and gravity tests, as well as tests of the tree and the error-scaling of the 
code. The iV-body component was tested with some 3-body examples that had known solutions. In this paper, we present 
tests of the combined scheme to demonstrate that the conservative equations of motion derived in Section [2.21 are correct, 
and that the scheme does not exhibit any unexpected numerical effects. We perform a simple test investigating the scattering 
between gas and star particles. Using star-gas Plummer spheres, we test the expected error scalings of the SPH and iV-body 
components and the combined scheme, as well as the stability of the star-gas Plummer spheres. Finally, we perform a simple 
test of colliding star-gas Plummer spheres. In all tests we use an adiabatic equation of state in which heating and cooling are 
only due to PdV work by expansion or contraction, thus enabling us to test the energy conservation of the code. We work in 
dimensionless units throughout, where G = 1. 



3.1 Star-gas particle scattering 

Our hybrid SPH/iV-body method enables investigation of stellar dynamics within a gas potential that may be time-evolving 
and irregular. For this goal, we must clearly understand the origin and impact of any numerical effects on the results of 
our simulations. One important difference between gas-only interactions and star-gas interactions in this scheme is that star 
particles are allowed to 'pass through' SPH gas particles, whereas artificial viscosity will prevent gas particle penetration by 
forming a shock. Even though the gravitational interactions between stars and gas are smoothed, the stars are still interacting 
with a gravitational field defined by discrete points and thus can be deflected by those points. Therefore, gas particles can in 
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principle scatter star particles significantly if the resolution is too coarse. This test is designed to investigate how significantly 
gas particles can scatter star particles and to help determine resolution criteria to prevent significant numerical scattering. 

For point particles obeying Newton's gravitational law, the scattering angle of a star of mass m a due to a hyperbolic 
encounter with a gas particle of mass m g in the centre-of-mass frame is given by 

A #d = 2mg tan-if^+^ggggi (31) 
m„ + m a \ bv 2 I bv 1 



(e.g. Binnev fc Tremaine 20081 ) where b is the impact parameter, v is the relative tangential velocity at infinity, and the 



approximation is for small deflections. For smoothed gravity, we would expect that the net deflection angle would be reduced 
by smoothing, and for the dependence on 6 to be fundamentally altered since the gravitational force for neighbouring particles 
tends to zero as the distance becomes zero. Fo r the M4-kernel, the gravitational force reaches a maximum at around |Ar| ~ 
0.8 h (see Figure 1 of Price fc Monaghianl 2007 ) and then decreases to zero (instead of cx 1/r 2 ). Therefore, an approximation 



to the maximum possible scattering angle can be obtained by setting b ~ h in Equation 13 llrl . Since the smoothing length is 
dependent solely on the gas particle positions, the strength of the star-gas interaction becomes solely a function of relative 
velocity, v. From Equation 1311 we can define a critical relative velocity, i) C rit, where the interaction results in a significant 
deflection angle, where 

= f^V /2 . (32) 



In this form, the deflection angle is simply A#d = v 2 lit /v 2 . Therefore, a simple resolution condition that ensures star-gas 
scattering is negligible, i.e. A6r> <C 1, is given by v 2 2> v 2 lit . This application of this criterion to various astrophysical 
scenarios will be discussed later in the paper. 

In order to test the validity of this assertion, we perform simulations of a single star that is moving with velocity v 
through a periodic gas cube of side-length L and uniform density p where the gas velocity is fixed to zero everywhere. If the 
gas density field is perfectly uniform (i.e. in the continuum limit where N g — > oo), then the net gravitational force due to the 
gas will be zero and the star will simply move with constant velocity without any deviation or deceleration. If the density 
field is not perfectly uniform, as is the case when represented by a discrete set of particles, then small deflections will alter 



the path of the star. In our test, the gas particles are first relaxed to a glass (See lHubber et al.ll201ll . for details) which is 



the most unifor m arrangement of par ticles used in typical simulations. We use periodic boundary conditions combined with 



Ewald gravity (|Hernouist et al.lll99lf ) to produce a net zero gravitat i onal f ield in the gas, with the exception of the small 



deviations due to the smoothed particle distribution. Fellhauer et al. ( 2000h performed a similar scattering test to quantify 



the effects of numerical relaxation, although in the context of large-scale galaxy simulations. 

The initial velocity v of the star particle is varied between simulations. We test v = jVcrit, ^Wcrit, v C rit, 2 w cr it , 4t»crit, 
8 Vcrit, and 16-Ucrit- For each case, we place the particle at a random position in the gas cube in order to remove the systematic 
effects of using the same glass arrangement. We simulate 6 realisations for each of the selected velocities and take the average 
and standard deviation of the results. Each star particle moves along its trajectory until t = L/v, i.e. the crossing time of 
the cube. This ensures, in so far as is possible, that during a simulation a star encounters approximately equal number of 
gas particles for all choices of initial velocity. To quantify the effect of the scattering on the star particle, we measure the 
deflection angle, A8u, i.e. the angle of the particle's trajectory relative to the initial velocity along the x-axis. The results of 
these tests can be seen in Fig. [2] 

For v > Wcrit, increasing the star particle velocity decreases the net effect of scattering due to the gas particles. In this 
regime, the scattering angle, A#d, falls as v~ 2 , the same dependency as suggested by Eon. 13 1 1 with b ~ h. We notice that 
the net average scattering angle is lower than expected (Figure [2] dashed line) due to the effects of smoothing. We note the 
net deflection is the accumulation of several (~ 10) deflections, not necessarily in the same direction hence it will 'random- 
walk' with each deflection. The error bars are due to the combination of the random-walk errors and the impact parameter 
dependence (which is reduced, but not eliminated). 

However, for v ^ «crit> extremely strong scattering occurs and dominates the dynamics of the star. At these velocities, the 
power-law relation between the scattering angle and velocity is broken since the interaction is now effectively a parabolic or 
elliptical interaction instead of a hyperbolic encounter. The final deflection angle is almost random due to the strong nature 
of the perturbations. The kinetic energy of the star is smaller than the gravitational potential energy, even accounting for 
smoothing, and therefore a star can in principle become trapped and bound to individual SPH particles. We note that this 
would not necessarily happen in a more realistic astronomical simulation because the gas particles in this test are static (so the 
star moves as a test particle), and therefore cannot be scattered off the star themselves. At such velocities, the coarseness of the 



2 A more rigorous derivation involving the form of the force kernel would reveal the exact dependence of the scattering angle on b and 
v. However, setting b = C h gives the same qualitative behaviour averaged over all impact parameters where C is some multiplicative 
constant of order unity. 
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Figure 2. Deflection angle A$rj versus star particle velocity (normalised by the critical velocity). As v increases above u cr i t , the star 
particle is increasingly less scattered. The scattering angle is roughly proportional to v~ 2 . However, for v ^ v cr a, there is a clear change 
in behaviour. Scattering is strong, and the scattering angle no longer follows the same power-law trend with particle velocity. 



SPH particle distribution clearly introduces potentially serious numerical effects which could corrupt any hybrid simulation. 
In future sections in this paper, we discuss possible gas resolution criteria as a means of avoiding unphysical scattering effects. 



3.2 Stability of star-gas Plummer spheres 



The Plummer sphere (|Plummeilll91ll ) is commonly used in stellar dynamics simulations as it is described by simple, analytic 
formula. For a Plummer sphere of mass, M, and characteristic Plummer radius, a, the density distribution, p(r), is given by 



p(r) 



3M 
47ta 3 



1 + 



-5/2 



and the 1-D velocity dispersion, u(r), is given by 

-1/2 



V) 



GM 
6 a 



1 + 



(33) 



(34) 



We simulate Plummer spheres that are purely stars, pur ely gas, or a m i xture of stars and ga^f]. For setting up the stellar 
component of our clusters, we use the method outlined bv lAarseth et alj (1 19741 ). 

For the gas distributions, the Plummer model corresponds to a n — 5 polytrope. A polytrope is a self-gravitating gas 
who se equation of state o beys the form P = A"p 1+ ™ and whose density structure is a solution of the Lane-Emden equation 
(See Chandrasekhai 19391 . for an in-depth description of polytropes). For the n = 5 polytrope, the radial density distribution 
is of the same form as Equation 1331 Instead of setting a velocity field complimentary to the density field to support against 
collapse, a polytrope is supported by a thermal pressure gradient. The thermal energy of the gas is related to the velocity 
dispersion of the gas by equating it to the sound speed and then converting to specific internal energy by u(r) — a 2 (r)/(y — 1) 
where 7 is the ratio of specific heats of the gas. We note that the gas itself in our simulation does not need to obey a polytropic 
equation of state, only that the thermal energy distribution of the gas is set-up to mimic the pressure distribution of the 
equilibrium polytrope and therefore remain in hydrostatic balance. The gas responds adiabatically and therefore can heat by 
contraction or cool by expansion as it set tles or is moved around by the potential of the stars. The initial positions are set-up 
using the method of lAarseth et al.l (1 19741 ), but the thermal energies are set using the above equation and the initial velocities 
are set to zero. 

We simulate the evolution of (a) a gas-only n = 5 polytrope, (b) a star-only Plummer sphere, and (c) a 50-50 mixture (by 
mass) of a star-Plummer sphere and a n = 5 gas polytrope. Gas-only simulations are conducted with N g = 5, 000 gas particles 
of total mass M g , and star-only simulations with N s — 500 star particles of total mass M s . For each case, M = M 3 + M g — 1 
and a = 1 using dimensionless units (where G = 1). Mixed star-gas simulations have either either N„ = 100 or 500 star 



3 Note that Plummer spheres are formally infinite in extent. In practice we truncate our Plummer spheres at a radius of 20 a. This 
means that they are not in exact equilibrium; however this has a negligible effect on their evolution 
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Figure 3. The fractional global energy error, AE/E, as a function of timestep multiplier, t for (a) SPH gas particles in a n = 5 polytrope 
integrated with a 2nd-order Leapfrog, (b) iV-body star particles in a Plummcr sphere integrated with a 4th-order Hermite, and (c) SPH 
gas and star particles in a Plummer/Polytropc combination where the gas is integrated with a 2nd-order Leapfrog and the stars with a 
4th-order Hermite, with and without the x correction term derived in Section I2.2I and Appendix [A] Also plotted for guidance are the 
2nd- and 4th-order scaling expected for the SPH, A^-body and hybrid simulations. 



particles and either 10 x or 100 x as many gas particles. In all cases we use equal-mass star particles (m a = M 3 /N 3 ) and equal- 
mass gas particles (m g = Af g /iV 9 jf]. The smoothing length of the stars in all simulations is h — 0.0001a. Each simulation is 
run for 40 crossing times, where we define the crossing time here to be i CR = a/a(r = 0) = ^pSaFjGM = 2.45 code time 
units. 



3.2.1 Resolution criteria 

Following the ideas discussed in Section \3. 11 we can determine the resolution requirements of an equilibrium Plummer sphere 
to significantly reduce the effects of unphysical star-gas scattering. Let us assume that N g gas particles account for a fraction 
/ of the total mass, i.e. M g = f M = N g m g where each gas particle has mass m g . The central gas density is po — 3 M g /A it a 3 , 
and the central velocity dispersion is a% = G M/6 a. We can then calculate the critical resolution velocity at the centre of the 
Plummer sphere as 



h J \ rja J \^4 7r 

where we have substituted Equation [T] for h and used the above expressions for po- 

In order to avoid catastrophic numerical scattering of star particles off gas particles, star particles must be moving at 
velocities significantly larger than the critical velocity, i.e. <tq 3> v c . This leads to the following resolution criterion for the 
total gas particle number, N, in the Plummer sphere as 



4 Equal-mass star particles is obviously a simplification for the purposes of our tests. However equal-mass gas particles should be used 
to reduce SPH noise and prevent particle clumping. 
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assuming the typical value of r\ = 1.2. For a Plummer sphere consisti ng of equal ga s and stellar mass (/ = 0.5) we find that 
N g 3> 5. For a gas-dominated Plummer sphere (/ = 1.0), N g ^> 15. iHubber et alj (|201lf ) demonstrated of order a hundred 
gas particles could only crudely represent an equilibrium polytrope, with reduced central density and a smaller radius. Such 
structures require of order thousands or tens of thousands of gas particles to adequately resolve the density structure of the 
polytrope. Therefore, it is clear that we require N > 1000 to resolve the density field, at which point the gas resolution is also 
sufficiently high to prevent serious star-gas scattering events from corrupting the simulation. 

We note that the stars are moving with a range of velocities, below and above the mean velocity dispersion, ao- No 
matter how high the resolution, there will always be a number of stars at some instant moving less than the critical velocity. 
Therefore, we cannot completely eliminate unphysical scattering in this case, but we can only reduce it to some acceptable 
level by using a reasonably high resolution. 



3.2.2 Energy conservation 

In order to test the energy conservation properties of hybrid code, we investigate how the fractional global energy error, 
AE/E, varies as a function of the timestep size. Instead of selecting a global, constant timestep, we allow an adaptive global 
timestep, and then vary the timestep multiplier, r (to which we equate all other timestep multipliers, 7 COUR , 7 ACCEL , 7 B nergy 
and 7 S , for this test), which determines the adaptive timestep size. Therefore we consider how the energy error scales with r, 
which should scale in the same way with a constant stepsize. 

Figure [3] shows the energy error scaling as a function of timestep multiplier for (a) gas-only polytrope, (b) star-only 
Plummer sphere, and (c) the star-gas mixture as a function of r. For the gas-only cluster, the SPH integration scheme is 
a 2nd-order Leapfrog kick-drift-kick. Therefore we expect the error to scale as |A_E| ~ 0(dt 2 ) ~ 0(t 2 ). In Figure EJa), we 
can see that the energy error (filled black circles) agrees very well with this expected scaling with only a small deviation 
from the added guideline (solid red line). For the star-only cluster, the JV-body integration scheme is a 4th-order Leapfrog 
scheme, and therefore we expect \AE\ ~ 0{dt A ) ~ 0(t 4 ). Figure [3jb) shows that we obtain similar scaling to this. In fact, 
we obtain slightly steeper scaling relative to the expected 4th-order. For the star-gas cluster, a combination of a 2nd-order 
scheme with a 4th-order scheme should in principle give either 2nd-order or 4th-order erros depending on whether gas-gas, 
star-gas or star-star interactions are dominating the error. Figure [3jc) shows that we get 4th-order scaling (filled black circles) 
which suggests that star integration scheme dominates the error in this simple case, despite the formally larger error of the 
gas integration scheme. We also plot the energy errors for they star-gas cluster with and without the new \ correction term 
(Equations 1 121 fe !15[l introduced in this paper. Without x, the equations of motion are non-conservative, and the energy error 
(of order 10~' j ) is dominated by this rather than integration error. 

As well as integration error, block timesteps and gravity tree errors are other major sources of error in this scheme, and 
will likely dominate over the integration error in most practical simulations depending on the parameters chosen. The error 
in the tree can be controll ed by an appropriate choice of mul tipole-acceptance cr iteria (MAC), such as the GADGET-style 



MAC (|Springel et al.l 12001! ) or the SEREN Eigenvalue MAC (|Hubber et al.ll201lh , both of which can place an upper-bound 



on the force error due to individual cells and hence control the net tree force error and indirectly the global energy error. 



3.2.3 Plummer sphere stability 



One of the aims of this hybrid scheme is to accurately follow the stellar dynamics of a system in a live gaseous background. 
The most significant difference between stars and gas in this context is that stars are collisional particles in the sense that 
they can be subject to strong two-body interactions, whilst gas is collisional in the sense that it can form sho cks. 

A low-iV stellar Plummer sphere can rapidly evolve due to two-body scattering. lAarseth et al.l (1 19741 ) showed, both 
numerically and through semi-analytical models, that the effect of two-body relaxation in a Plummer sphere is to redistribute 
energy ejecting some stars and causing the contraction (eventually to core collapse) of the remaining cluster. lAarseth et al. 
l| 19741 ) traced this with the Lagrangian radii of the stellar clusters showing the contraction of the inner Lagrangian radii and 
the expansion of the outer Lagrangian radii due to energy conservation as stars are ejected. They showed that the 50 per 
cent Lagrangian radius stays relatively constant throughout the evolution. The Lagrangian radii evolve significantly on the 



two-body relaxation timescale, i ro i a 
0.1 N s 



: , given in terms of the crossing time, i c 



^rcla 



In N a 



■t c 



(37) 



where N s is the number of stars (See Binnev fc Tremaine 200Sh 



The evolution of the star-o nly Plummer sphere (Figure|4ja)), containing 500 stars, shows the same qualitative behaviour as 
found bv lAarseth et al.l (|1974l ) demonstrating that our pure A r -body integration scheme is correctly capturing the qualitative 
effect of 2-body encounters over the expected timescale (for N g — 500, the relaxation timescale is t ro iax ~ 10£ cr oss using 
Eon. I37[) . As observed by Aarseth (1974) and Aarseth et al. (1974), the 10 per cent Lagrangian radius shrinks (towards core 
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Figure 4. The time evolution of the 10%, 50% and 90% Lagrangian radii for Plummer spheres of M = 1, a = 1 (in dimensionless units) 
containing (a) 500 stars only, (b) 5, 000 gas particles only, (c) 5, 000 gas particles and 500 stars, (d) 50, 000 gas particles and 500 stars, 
(e) 1,000 gas particles and 100 stars, and (f) 10,000 gas particles and 100 stars. 

collapse), the 90 per cent Lagrangian radii expands (due to ejections), and the 50 per cent Lagrangian radii stays roughly 
constant. 

The gas-only Plummer sphere (Figure HJb)), containing 5,000 gas particles, is observed to evolve slightly over the first 
few crossing times as it settles into equilibrium but soon the Lagrangian radii stay almost constant over time. 

Of far greater interest is the behaviour of a mixed star-gas Plummer sphere. We run two sets of simulations: 100 star 
particles with either 1000 or 10000 gas particles, and 500 star particles with either 5000 or 500000 gas particles. We use two 
different gas resolutions for each case to help determine if the simulations are converged. The evolution of the 10, 50 and 90 
per cent Lagrangian radii of both the stars and gas in each case are shown in Figure [4] 

We notice, for both N s = 500 and N s = 100, that the Lagrangian radii for the stellar and gas components evolve in 



A hybrid SPH/N-body method for star cluster simulations 13 




X 



Figure 5. Supersonic collisions between gas-dominated (90 per cent gas by mass) Plummer spheres at (a) the initial state; (b) just after 
the collision; (c) the end of the simulation. Each Plummer sphere has 5000 gas particles and 200 equal-mass star particles. The initial 
crossing time of a Plummer sphere is 2.45 code units and the time is measured in code units. Stars are shown by white dots, the colour 
table shows the column density of the gas in code units. 

opposite directions. The stellar component is altered somewhat from the star-only case where it appears to shrink for the 10 
per cent and 50 per cent Lagrangian radii (with the 90 per cent remaining fairly constant). Conversely, the gas appears to 
expand at all radii, and on a timescale comparable to the stellar 2-body relaxation timescale. Therefore, there is a transfer of 
energy from the stars to the gas, allowing the gas to heat and expand, and conversely the stars lose energy and contract. 

Most importantly for this paper, the results converge for different gas particle numbers (with the same number of star 
particles). For N s = 500, we use both 5,000 (Figure HJc)) and 50,000 (Figure[4jd)) gas particles. The evolution of both the 
N s = 500 simulations is basically identical. There is some deviation between the N s = 100 results at late times in different 
N g backgrounds. This is due to low-iV s noise and the slightly earlier 'core collapse' of the N g = 10, 000 simulation. 

The reader might notice that the behaviour of the stars in the N s = 500 star-only simulation is somewhat different to that 
of the stars in the N s = 500 simulations with gas. This is an interesting physical (not numerical) effect due to the presence 
of gas. We will return to the physics and astrophysical implications of this behaviour in the next paper. For now, however, 
we will simply use these simulations to illustrate the convergence of the results for different numbers of gas particles but the 
same number of star particles. 



3.3 Star-Gas cluster collisions 

As a simple qualitative test of the hybrid code's ability to model more complex star-gas systems, we perform a small suite of 
simulations of the head-on impact between two star-gas Plummer spheres. We create two star-gas Plummer spheres following 
the procedure described in Section 13.21 We collide the Plummer spheres at a velocity u co ii, such that the collision occurs 
either subsonically or supersonically for all gas particles. We also consider Plummer spheres that are gas-dominated, and 
star-dominated. We therefore expect strong differences in the behaviour of the gas and stellar dynamics between the subsonic 
and supersonic tests, and also between the star and gas-domainted cases. 

Each Plummer sphere contains N g = 5, 000 equal-mass gas particles and N a — 200 equal-mass star particles and is set-up 
in the same way as in Section [3.21 For gas-dominated cases, 90 per cent of the mass is in gas and for star-dominated cases, 
90 per cent of the mass is in stars (therefore the relative masses of star and gas particles are different by a factor of ~ 100 
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Figure 6. Supersonic collisions between star-dominated (90 per cent stars by mass) Plummcr spheres. Otherwise the same as fig. [5] 



between the two cases). For subsonic collisions, the relative collision velocity is 0.5 (in dimensionless units), and for supersonic 
collisions, the relative collision velocity is 10. 

Figures [5] and [6] show the results of supersonic collisions of a gas-dominated and a star-dominated collision respectively. 
In the gas-dominated supersonic collision (Fig. [5J, the gas forms a shock around x = where the gas is heated up and is 
compressed to higher densities. Meanwhile the stars pass through the shock front and also through the stellar-component of 
the other cluster. Since the relative velocity of the two clusters is much greater than the individual velocity dispersions, the 
effects of two-body encounters are negligible and the two clusters pass through each other almost unperturbed. As the gaseous 
and stellar components have decoupled in the collision, the gas-free stellar clusters are now unbound as they have had 90 per 
cent of their initial mass (the gas) removed. Therefore they expand as their velocity dispersion is too high to maintain their 
initial dense configuration ; the cluster eventually dissolves over several crossing times (this is analogous to gas expulsion, see 
Goodwin fc Bastian 20061 ) . In the star-dominated supersonic collision (Fig. [6]), the gas again shocks and decouples from the 



stars. However, as the stars dominate the potential, the cluster is only slightly super-virial (| < Q < 1) and the degree of 
expansion whilst readjusting to the new potential is small. Therefore, the two stellar clusters continue with almost no effect 
from the collision and the stripping of their gas. 

In Figures and [S] we show the results of subsonic collisions of gas-dominated and star-dominated clusters respectively. In 
both cases the two clusters merge as the gas components merge and the stellar components can respond to the interaction as 
they interact on a timescale comparable to their crossing times. The shapes and details of the resulting clusters are different; 
the star-dominated merger showing a more elongated final appearance (compare the last panels of Figures [7| and [8]) . This is 
explained as the post-shock velocity anisotropy of the stellar velocity field dominates in the star-dominated case, whilst in the 
gas-dominated merger the potential is dominated by the gas allowing significant violent relaxation of the stellar component 
in the spherical gas potential. 

The behaviour of the supersonic and subsonic collisions in these simulations is physically reasonable. We also track energy 
conservation and find it is very well conserved ( AE/E = 1.6 x 10 -6 — 1.4 x 10 -5 for the four simulations) . As with the static 
Plummer test, we will explore the physics of star-gas cluster collisions in more detail in a future paper. 
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Figure 7. Subsonic collisions between gas-dominated (90 per cent gas by mass) Plummer spheres. Otherwise the same as fig. [5] 
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Figure 8. Subsonic collisions between star-dominated (90 per cent stars by mass) Plummer spheres. Otherwise the same as fig. O 
4 DISCUSSION 

Our main motivation for developing this new hybrid method is to provide an intermediate step between detailed hydrodynam- 
ical simulations and pure iV-body simulations. In particular, we wish to investigate the dynamics of stars in gas on > 0.1 pc 
to pc-scales in GMCs and young clusters. The main advantage of our method is that the dynamics of stars within a live gas 
background can be followed at relatively low computational expense compared to full hydrodynamical simulations (hours or 
days on desktop computers compared to months on HPC facilities). We discuss a number of important practical considerations 
for preparing simulations using the hybrid code, as well as planned future developments and uses of the code. 



4.1 Suppressing fragmentation 

We reiterate that this method is not intended to act as a replacement for full hydrodynamical simulations which aim to 
model the star formation process itself, i.e. the fragmentation of molecular clouds into prestellar cores and finally multiple 
protostellar systems. The fragmentation of gas into stars is a complex hydrodynamical problem, involving much additional 
physics such as radiation transport, chemistry, and (non- ideal) MHD. We suggest that fragmentation should be artificially 
supressed for a number of reasons. Firstly, we wish to avoid the complex physics and computational expense related to full 
hydrodynamical simulations since we are only currently interested in the global effects of the gaseous gravitational potential 
on the iV-body evolution of the cluster. Secondly, if fragmentation occurs then sink particles should be introduced and it is 
unclear how to mix sink particles with iV-body star particles (one will be softened and interacting hydrodynamically with the 
gas, the other will not despite both representing a star). 

Therefore we suggest that the resolution be kept as low as possible and that the equation of state be designed to supress 
star formation and keep densities low. We note that this also has the advantage of avoiding the formation of discs around stars 
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which are again complex, computationally expensive objects to model. We suggest that hybrid simulations are designed to 
avoid the regimes in which these processes are important, ie. sub-core (< 0.1 pc) s cales and high gas densi ties (> 1CP 13 g cm~ 3 ). 
The normal resolution criteria for SPH simulations of star formation is the lBate fc Burkertl (| 19971 ) criter ia in which the 



(min imum) Jeans mass must be resolved i n order to cap t ure g ravitational fragment ation (or in AMR, the iTruelove et al 



19971 . Jeans length criterion). As shown by iHubber et al.l (|2006l ). failure to meet the lBate fc Burkertl (|1997l ) criteria means 



that fragmentation is suppressed - this is important as it means that low SPH resolution results in no fragmentation rather 
than artificial fragmentation. 

We suggest tha t gas be kept at densities low er than 1CT 



•10" 



g cm , well above the critical density for fragmentation 



(i.e. 10 13 g cm 3 , Masunaga fc Inutsuka 2000l ) . Such densities are in the roughly isothermal regime, and so have a simple 



equation of state (e.g. Jappsen et al.l 2005 ). We suggest that the equation of state be modified (say by artificial heating) to 
keep densities low. The potential to form cores is probably desirable, but not to follow their collapse and fragmentation. 



4.2 Avoiding unphysical scattering 

As discussed and demonstrated in the star-gas scattering test (Section 13. ip . a sufficiently high gas resolution is required to 
avoid the unphysical scattering of stars by gas particles. For a general situation, where there is a group of stars moving through 
a cloud of gas, but no equilbrium has been established, then the velocities are not linked in any way to the density of the 
gas and it is difficult to establish a simple criterion for the required gas resolution. There are various scenarios where we can 
establish a link between the star velocity and the gas density. For example, equilibrium Plummer spheres, we were able to 
derive a resolution condition for equilibrium Plummer spheres on the number of gas particles which was only a function of 
the gas mass fraction, / (Eqn. l36]) . We note that this condition is only strictly true for equilibrium clusters where the velocity 
(or velocity dispersion) of the stars is well-known. For non-equilbrium simulations, we may need to use further information 
infered from the initial conditions, such as the initial virial ratio, to infer how the velocity relates to the density, and hence 
to the required resolution of the gas. 

One other special scenario is a gaseous cluster with emedded primoridal binary and multiple systems which may be 
unphysically disrupted due to star-gas scattering. Consider the simple case of a binary star containing two stars of mass 
mi and ui2 in a circular orbit of separation a BIN , with an orbital velocity w BIN = \/G(m\ + m2)/a BIN . To avoid artificial 
scattering from potential destroying the binary, we require that w BIN 2> w C rit • Assume that a binary is located at the centre 
of a pure gas plummer sphere, consisting of N gas particles, of total mass M, and scale-length a. The critical velocity of the 
central gas is given by Equation [35] Rearranging leads to the resolution condition, 

3/2 

(38) 



N > 



2/ x 3/2 / 3 X l/2 



r\ ) \4-7r 



M 



(mi + m,2) 



Inserting in reasonable values of a = 0.1 pc, a BIN = 10 AU, M = 10 Mq, and mi + 7712 = 0.1 Mq, we find N 3> 200. Therefore, 
for a smooth distribution of gas, > 10, 000 is sufficient to avoid star-gas particle scattering. We note that wide binaries 
would be most sensitive to star-gas scattering in comparison to tight binaries. Therefore, using wider-separation primordial 
binaries would require higher resolution. However this scenario is highly idealised compared to more practical scenarios. A 
more irregular gas distribution may have higher values of u cr it in high density pockets resulting in stronger scattering and 
therefore, more stringent resolution requirements. 



4.3 Future applications and code developments 



There are many possible future astrophysical applications for our hybrid method. Our first follow-up papers will address 
stellar dynamics in small- A'' star-gas groups, and the collisions and mergers of such objects as touched upon in this paper. 

In the longer term, we plan to simulate larger more complex systems (dynamics in turbulent gas or fractal distributions) 
both as simple numerical experiments and to compare with observations. We plan to add simple stellar heating and gas 
cooling prescriptions in order to advance on the simple adiabatic EOS we have employed here. For feedback and gas dispersal 
simulations, we will add some simple feedback formulations, like gas-heating from sup ernovae, as well as mechanical winds 
and UV radiation using the HEALPix-based algorithm already implemented in SEREN jBisbas et al.ll2009l ). Simple accretion 
models can be added to allow stars to accrete from the gas. 

One important algorithmic addition we are currently implementing in the code is a more sophisticated N-body inte- 
grator that will allow efficient evaluation of computationally expensive sub-sytems such as tight binaries and 3- or 4-body 
encou nters. We are implementing an adaptive n earest-neighbour tree, similar to that used in the STARLAB N-body code 
suite ( Porteeies Zwart et al. 2001) and MYRIAD ( Konstantinidis fc Kokkotas 2010l ). to decompose the stars into sub-systems 
and then use a higher-order integrator, such as the 6th or 8th order Hermite integrators, to more accurately integrate the 
sub-system. This addition will allow us to model clusters containing primordial binaries and higher-order multiple systems, or 
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clusters that form hard binaries. Details and tests of any additional physics and optimisations will be explained in subsequent 
papers that introduce them. 



5 CONCLUSIONS 



We have presented a new hybrid SPH/iV-body method within the SPH code SEREN (jHubber et al.ll2011f ). Using conservative 
SPH and 4th-order iV-body integrators, this scheme conserves energy extremely well with an adiabatic equation of state. We 
have presented a number of tests of the code showing that it works as expected in a number of simple situations. We will 
use this code in future to explore stellar dynamics in a live gas background to investigate problems involving star cluster 
formation and evolution. 
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APPENDIX A: DERIVATION OF CONSERVATIVE SPH EQUATIONS WITH STARS 



Following iPrice fc Monaghan! 1)20071 ). we derive the SPH equations of motion, for both the gas and star particles, using 
Lagrangian mechanics. For a set of N g gas particles with labels b = 1, 2, N g and N a star particles with labels i = 1, 2, N s , 
then the Lagrangian becomes 



C = \ Yl ^ V * + \ mtV i + £ GRAV ( A1 ) 

6=1 i = l 

where £ GRAV is the gravitational contribution to the Lagrangian given by 



g N g N g g n b N B N g N s 

^grav = Y2 m bm c 4> bc (h bc ) - — 22 Y2 m i m j ( l>ij(hij) - G^2 Y2 m b m i4>bt(hbi) (A2) 

6=1 c=l i=l j = l 6-1 i=l 

where h bc = \ (hb + h c ) is the mean smoothing length of particles b and c. The equations of motion can be obtained by 
inserting the Lagrangian into the Euler-Lagrange equations, 

- » 

The Lagrangian is symmetric in terms of interaction terms between the gas and star particles. The only difference lies in the 
method of calculating the smoothing length which leads to different forms of the equation of motion for both stars and gas. 



Al Gas particles 

First, we derive the equation of motion for a general gas particle labelled a by taking the derivative with respect to its position, 
r 0) i.e. 

\' iT \ - is c \ . f dcf>bc (dh b dp b dh c dp c \ 

Pbc{nbc)T-bc(0ba-O C a) + -^=- Tj T, h Tj T : 

2 dh bc V "Pb or a dp c dr a J 



Ng Ns 

— mb rrii 

6=1 i=l 



^ - -f EE m » m » 



4>bi(hbi) Hi Sba + 



1 d(/)bi dh b dpb 

2 dhu dpb dr a 



(A4) 
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We note there is no contribution from the star-only term in the Lagrangian since there is no dependence on the position of 
any gas particles, i.e. r a - Substituting the expression for dp/dr, i.e. 



dp b 1 dW bd (h b ) 

fr- a = TT b 2^ md —d^ {Sba - Sda) - 



(A5) 



where Q b is given by Equation [13] we obtain 



dC r 



G 



Ng N. 



mb mc ^'bdhbc) he (Sba - S ca ) - m; 4>'bs{hbs) hs S ba 



6=1 c=l 

Ng Ng Ng 



6=1 c=l d=l 

Ng N s Ng 



6=1 i=l 

1 dh b dWbd(h b 



G \ -> \ -> \ -> Ocpbc 

~^^^ mbmcmd ^\n b d Ph dr a 



(06a — Oda) + 77- ~ (f>ca — Oda J 



G v-^ vA ^ depbs 1 dh b dWbd(h b ) , . . 



6=1 i=l d=l 



(A6) 



Expanding out the Kronecker delta functions and simplifying, 



dC r 



dr a 



N„ N„ jv, 

G G 

— ^2,m a m c <j)' ac (h ac ) f oc + — ^2 m b m a (f>' ba (hb a ) r ba - G^2,m a mi <j)' ai (h a3 ) r ai 



c=l 

g «j 



JV„ N 

c=l d=l 

tEE 1 

6=1 d=l 

N, Ng 



G 



rrii rrid 



d<j> ac 1 dh a dW ad (h a ) , G 1 9/i dW ba {hb) 

+ T / y rn b m c m a -=— 77- ■- 

d/i 6c n 6 9p6 dr a 

d(j) bc 1 <9fc c dW ea (ftc) 
dh bc fi c 9p c <9r a 

<906» 1 dh b dWba{h b ) 
dh u fi 6 dp b dr a 



dhac 




dv a 


906a 


1 dha 


dWad(ha) 


9/l6a 


fia <9p a 


dr a 




1 a/i a 


dWad(h a ) 


dh ai 


fia 9pa 


dr a 




N a 





6=1 c=l 

G ^ ^ 

+ jEE™'™ 1 ™ 

6=1 c=l 



+ 1EE 



-iEE m « 

i=l a=l 

— G^ma m 4>' a b{hab)r a b - G 2^ m a m, (fi' as (h a i) r a % 



mb mi m a 



G 



6-1 



Ca 0W„ 6 (ft a ) , Cb OTM/*) 



+ 



fi a <9r a n 6 <9r, 



G 



'"6 



Xa gWai(/ta) X6 gWqfr(M 

fi a dv a Q b dv a 



(A7) 



where £ a (cf. iPrice fc Monaghanll2007f ) and Xa are defined by 



J dh a v~> ^0a6 

9pa <%a6 



(A8) 



9/la <90ai 

9 P<* ~i dhai 



(A9) 



Substituting into the Euler-Lagrange Equation (Equation I A3|) . we obtain the equation of motion for SPH gas particles, 



a a = ~G (j>' ab (hab) i ab -G^ j uij <j>' a3 (h aa ) r aa - --■ 



'"6 



(Ca + Xa) aWa6(fta) . (Cb + X6 ) 9 W a 6 (ft 



fia 9r a 



+ 



S>6 



9r a 



(A10) 
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A2 Star particles 

Similarly for star particles, we derive the equation of motion for a general star particle labelled a by taking the derivative of 
the gravitational component of the Lagrangian with respect to its position, r a , i.e. 

dC G Ns Nb N " Nb 

— g y RAV = --2^2^2 mam i [<t>ij(hij)rij (S ia - S ja )] - G^2^2m b mi y>' bi (hu) ? u (-&<»)] 

a i=l j=l 6=1 i=l 

G N ' , - G N " - Ns 

= — n yi"»a rrij fiatjhgj) r aj + — ^ m, m a (j}' aa {h ia ) i ia + G ^ m b m a (j}' ba (h ba ) f ba 

3=1 i=l 6=1 

= -G}^m a mict>' ai (hai)r ai — G}^ m a m b <t>'ab{h a b) r a b (AH) 

i=l 6=1 

Due to the stars having constant smoothing length, we obtain somewhat simpler equations than for the case of gas particles. 
Substituting into the Euler-Lagrange equations and renaming some summations for clarity, we obtain the following expression 
for the acceleration of star s, 

N s Ng 

a s = ~G^2 irn <f>' si (hsi) r 3i - m b <f>' sb {h sb ) r sb (A12) 



